2.1.1.1 Population numbers

These data will be provided by the synthetic population or derived from other data. In the code presented here, we simulate 5-year age and sex cohorts. Data for population may be in five-year age groups or one-year. For the example for Greater London, we derive 5-year age groups from GBD IHME data, however, we also provide a script if five-year age groups are to be derived from one-year age groups data. I left potential data sources below as a reference.

Data source: (1) National census; (2) Worldwide population and mortality data: http://www.mortality.org/ (mostly high income countries; and (3) Calculate from the Global Burden of Disease by the Institute of Health Metrics and Evaluation (GBD IHME) data (rates and numbers available from (http://ghdx.healthdata.org/gbd-results-tool).

2.1.1.2 Mortality rates

Mortality rates are needed per single year and sex. These data are available from GBD IHME, however, in age groups (1-4, 5-9, etc). We provide an interpolation script to derive in between ages rates (cubic spline).

Note that we need data for population numbers and all cause mortality rates for: (1) PMSLT and (2) Dismod II collection (more in Dismod II section). Population data from the synthetic population is used for the PMSLT (if available). For Dismod II, population and mortality data should be from the same source (GBD IHME).

2.1.1.3 Total years lived with disability rates per single year and sex.

These data is available from the GBD (http://ghdx.healthdata.org/gbd-results-tool) per 5-year age groups. We can use interpolation to derive between ages rates (script provided).

2.1.1.4 Disability weights (quality of life weights)

Disability weights (DW) are derived from disease specific years lived with disability (YLD) and disease specific prevalence by age group (5 years) and sex. Data for YLDs prevalence is obtained from the online GBD IHME data tool (http://ghdx.healthdata.org/gbd-results-tool). An age and sex specific-correction was introduced to counteract the effects of accumulating comorbid illnesses in the older age groups (Equation 1).

\[\begin{equation} \label{DW adjusted for total YLDs-YLDdPd1YLDt} (YLDd/Pd)/(1-YLDt) = DW adjusted for total YLDs \end{equation}\]

Where YLDd is the YLD mean number per age and sex for a given disease, Pd is the prevalence (as reported in GBD ) for a given disease by age and sex and YLDt is total YLD rate per age and sex.

2.1.2.1 Incidence and case fatality

For each of the modeled diseases the PMSLT needs incidence and case fatality rates per sex and one-year intervals. Data from the GBD IHME studies with Dismod II (free at https://www.epigear.com/index_files/dismod_ii.html) is used to derive internally consistent data and generate missing data. For example, the GBD studies provide data for incidence, prevalence and disease mortality, however, not case fatality. Other national level sources may also be explored/used, and compare with estimates produce from GBD data and Dismod II.

Dismod II inputs are: (1) population numbers and mortality rates and (2) disease specific inputs.

Population and mortality

Within Dismod II, each setting (e.g. country) has a collection that consists of population numbers (preferably the same as used in GBD IHME studies, due to the mortality envelop) and all- cause mortality rates (numbers and calculate rates). The GBD provides 5-year age groups that are acceptable input parameters for Dismod II.

Disease inputs by age group and sex

Each setting collection has a given number of diseases. Dismod II works with at least three of: case fatality, prevalence, incidence, mortality (disease), case fatality, remission, duration and the relative risk for mortality. So far, we have been assuming that remission is zero for chronic diseases, that is to say, when people become diseased, they do not recover. Special care should be taken with this assumption, as the GBD data assumes remission for some diseases, for example cancers, where after 10 years cases recover, except for long term sequelae. Since GBD now provides prevalence, incidence and mortality, it may be best to use all three as Dismod II input parameters to compare the effect of the remission assumption by the GBD for some diseases.

2.2.2.1 Global Burden of Disease data

First, we explain how to obtain the data, second additional processing to derive data not reported (population) and one-year age groups (original data is in five-year age groups) and last procedure to use Dismod II. Data from the Global Burden of Disease data in Figure 3 can be download from here: http://ghdx.healthdata.org/gbd-results-tool. Figure 4 is a screenshot of the GHDx.

Table 2 specifies the selections to do for each of the tabs in Figure 2.

Once the selections described in Table 2 are made, the option *Download CVS** in the GHDx website is selected. A prompt comes up asking for an email address. The data is sent to the designated email address (within minutes) in ZIP format, unzip and use the code below to read the data (4). Here, we selected data for Greater London and England. The aim is to compare and derive scaling factors as for most cities the data is not available from the GBD and country level data may be used and scaled to the city level. Note that the data input requirement for the PMSLT, except population numbers, is in rates. Therefore the scaling is to better reflect the burden of an area, this is a different issue than working with numbers (e.g. total mortality numbers, total YLDs numbers) as in the ITHIM approach.

4. Read GBD data

GBDdata <- read.csv("data/UK/gbd2016.csv", stringsAsFactors = F)

The following codes serves to sort out the GBD data to the inputs required for the life table, disease life table and Dismod II

These data should be used to generate the general life table and disease life tables (Figure 3).

5. Change all upper cases to lower case and delete () from variables.

GBDdata <- mutate_all(GBDdata, funs(tolower))
GBDdata$measure[GBDdata$measure=="ylds (years lived with disability)"] <- "ylds"

6. Create age categories index in GBDdata (mid age, to match cohort running), total of 20 age groups. These are the age cohorts to simulate.

# 22, 27, 32, 37, 42, 47, 52, 57, 62, 67, 72, 77, 82, 87, 92, 97
GBDdata$age_cat [GBDdata$age =="under 5"] <- 2
GBDdata$age_cat [GBDdata$age =="5 to 9"] <- 7
GBDdata$age_cat [GBDdata$age =="10 to 14"] <- 12
GBDdata$age_cat [GBDdata$age =="15 to 19"] <- 17
GBDdata$age_cat [GBDdata$age =="20 to 24"] <- 22
GBDdata$age_cat [GBDdata$age =="25 to 29"] <- 27
GBDdata$age_cat [GBDdata$age =="30 to 34"] <- 32
GBDdata$age_cat [GBDdata$age =="35 to 39"] <- 37
GBDdata$age_cat [GBDdata$age =="40 to 44"] <- 42
GBDdata$age_cat [GBDdata$age =="45 to 49"] <- 47
GBDdata$age_cat [GBDdata$age =="50 to 54"] <- 52
GBDdata$age_cat [GBDdata$age =="55 to 59"] <- 57
GBDdata$age_cat [GBDdata$age =="60 to 64"] <- 62
GBDdata$age_cat [GBDdata$age =="65 to 69"] <- 67
GBDdata$age_cat [GBDdata$age =="70 to 74"] <- 72
GBDdata$age_cat [GBDdata$age =="75 to 79"] <- 77
GBDdata$age_cat [GBDdata$age =="80 to 84"] <- 82
GBDdata$age_cat [GBDdata$age =="85 to 89"] <- 87
GBDdata$age_cat [GBDdata$age =="90 to 94"] <- 92
GBDdata$age_cat [GBDdata$age =="95 plus"] <- 97

7. Create age and sex categories to obtain population numbers. Population numbers from GBD are used in Dismod II. For the Life table (Figure 3), the numbers may be from the synthetic population. For now, the Life table is set up with population numbers derived from the GBD data.

GBDdata$sex_age_cat <- paste(GBDdata$sex,GBDdata$age_cat, sep = "_"  )

8. Convert string variables to numeric to do calculations.

GBDdata$val <- as.numeric(as.character(GBDdata$val))

9. Generate population numbers for Greater London in a new data frame (“GBD_population”). Note that there is data for England as well, which is used in a separate rmarkdown document (GBDCompare). Population numbers are derived from rates per 100,000 and total numbers of cases.

GBD_population <- filter(GBDdata, measure == "deaths", cause == "all causes", 
metric == "rate" | metric == "number" ) %>% select(metric, age_cat, sex, val, 
sex_age_cat, location)

10. Generate population numbers from given number of cases and rates per 100,000 people.

 for (i in 1:nrow(GBD_population)) {
    if (GBD_population$metric[i] == "number") {
      GBD_population$five_year_population[i] <- GBD_population$val[i] * 100000/
      GBD_population$val[i + 2]}     
    else {GBD_population$five_year_population[i] <- NA}
 }

11. Remove rows with zero

GBD_population <- GBD_population[!is.na(GBD_population$five_year_population),]

12. Keep relevant variables

GBD_population <- filter(GBD_population) %>% 
select(sex_age_cat, sex, age_cat, five_year_population, location)

13. Create data frames for Greater London to be later used for: a) interpolation of rates and b) PMLT cohorts.

GBD_population_GL <-  filter(GBD_population, location == "greater london") %>% 
select(sex_age_cat, age_cat, sex, five_year_population, location, sex_age_cat)

12. Check population total numbers

GreaterLondon <- sum(GBD_population_GL$five_year_population)

13. Generate data frames for Greater London with per person rates (per 100,000 in original data).

GBDGL <- filter(GBDdata, location == "greater london" & metric == "rate") %>% 
select(measure, location, sex, age, metric, cause, val, age_cat, sex_age_cat)
GBDGL$one_rate <- GBDGL$val/100000

14. For the life table, we need to mortality and total yld rates in one year age intervals. The original data is in five years. Thus, the following code is used to interpolate a single-year age distribution form a five-yearly distribution.

## Loop to generate interpolated rates for all cause mortality and ylds for males and females
## UPDATE WITH YLDS RATES ALL CAUSES ADJUSTED FOR ALL OTHER MODELLED DISEASES.

i_sex <- c("male", "female")
i_measure <- c("deaths", "ylds") #" (years lived with disability)")

for(sex_index in i_sex) {
  for (measure_index in i_measure) {
   
data <- filter(GBDGL, measure == measure_index, sex == sex_index, 
cause == "all causes") %>% select(measure, location, sex, age, metric, 
cause, val, age_cat, one_rate) 
assign(paste(sex_index, measure_index, "interpolated_data", sep = "_"), data)
x=data$age_cat
y=log(data$one_rate)


interpolation_func <- splinefun(x, y, method = "natural", ties = mean)

interpolated <- as.data.frame(interpolation_func(seq(0, 100, 1)))
age <- seq(0,100,by=1)
interpolated <- cbind(interpolated, age)
interpolated[,1] <- exp(interpolated[,1])
colnames(interpolated)[1] <- paste(measure_index)
## Add column with sex to create age_sex category to then merge with input_life table
interpolated$sex <- paste(sex_index)
interpolated$sex_age_cat <- paste(interpolated$sex, interpolated$age, sep = "_")
## Change name of column death to mx and ylds to pyld_rate to then merge 
## with input_life table

if (colnames(interpolated)[1] == "deaths")
    colnames(interpolated)[1] <- paste("mx")
  else
    colnames(interpolated)[1] <- paste("pyld_rate")

# Name data frame
assign(paste(sex_index, measure_index, "interpolated", sep = "_"), interpolated)
# plot(interpolated$age,interpolated$rate)

##Do graph with another layer for original rates in age groups for comparison purposes.
p <- ggplot(data = interpolated,mapping = aes(age, interpolated[,1])) +
              geom_line(aes(color = "Interpolated")) +
              geom_point(
              data = data, 
              mapping = aes(age_cat, one_rate, color = "Original")) +
              labs(colour="",x="Age",y="Rates", sep = " ") +
              labs (title = paste("Rates", sex_index, "all cause", 
              measure_index, sep = " "), size=14) +
              theme_classic() +
              theme (plot.title = element_text(hjust = 0.5))
              print(p)
  }
}

15. The following code brings together all the required data to create the inputs life table (Figure 3).

## Create empty data frame
input_life_table <- as.data.frame(matrix(0, ncol=2, nrow = 202))
## Give variables names (mx=total mortality rate, 
## yldsx = total years lived with disability rates )
names(input_life_table) <- c("age", "sex")
## Populate life table: Age, 0 to 100
input_life_table[1:101, 1] <- c(0:100)
input_life_table[102:202, 1] <- c(0:100)
## Populate life table: sex, male, female
input_life_table[1:101, 2] <- "male"
input_life_table[102:202, 2] <- "female"
## Create variable age_sex to match with population data. 
input_life_table$sex_age_cat <- paste(input_life_table$sex,input_life_table$age,
                                sep = "_"  )
##Populate life table: five_year_pop. 

input_life_table <- merge(input_life_table, select(GBD_population_GL, c(sex_age_cat, five_year_population)), by = 'sex_age_cat', all = TRUE) 
##Populate life table: mortality rates (from interpolated rates)
## Males mortality rates. NEED MATCHING NAMES
input_life_table <- merge(input_life_table, select(male_deaths_interpolated , c(sex_age_cat, mx)), by = 'sex_age_cat', all.x = TRUE)
input_life_table <- merge(input_life_table, select(female_deaths_interpolated ,c(sex_age_cat, mx)), by = 'sex_age_cat', all.x = TRUE)
input_life_table <- merge(input_life_table, select(male_ylds_interpolated ,c(sex_age_cat, pyld_rate)) , by = 'sex_age_cat', all.x = TRUE) 
input_life_table <- merge(input_life_table, select(female_ylds_interpolated ,c(sex_age_cat, pyld_rate)), by = 'sex_age_cat', all.x = TRUE)

## DISCUSS WITH ALI TO IMPROVE

input_life_table$mx <-  ifelse(input_life_table$sex == "male",
input_life_table$mx.x, input_life_table$mx.y)
input_life_table$pyld_rate <-  ifelse(input_life_table$sex == "male",
input_life_table$pyld_rate.x, input_life_table$pyld_rate.y)
## Drop redundant variables
input_life_table <- subset(input_life_table, select = -c(mx.y, mx.x, pyld_rate.x, pyld_rate.y))

## Cross check population numbers
sum(input_life_table$five_year_population, na.rm = TRUE)

## [1] 8778396

## Sort data by sex and age to use in life table function. 
input_life_table<-input_life_table[order(input_life_table$sex, input_life_table$age),] 

16. Generate baseline life tables from input_data_frame using run_life_table function. Life tables are generated for each age (5-years) and sex cohort. It is assumed that current observed rates of mortality and disability are going to be observed in the future.

i_age_cohort <- c(22, 27, 32, 37, 42, 47, 52, 57, 62, 67, 72, 77, 82, 87, 92, 97)

general_life_table_list_bl <- list()

index <- 1

for (age in i_age_cohort){
  for (sex in i_sex){
    # cat("age ", age, " and sex ", sex, "\n") #Uncomment to see index
    general_life_table_list_bl[[index]] <- run_life_table(in_idata = input_life_table,
    in_sex = sex, in_mid_age = age)
    index <- index + 1
  }
}

## Uncommnet to check life table list
# View(general_life_table_list_bl[[2]])

17. Generate baseline disease life tables using run_disase function. Inputs of the disease life tables are: incidence, case fatality (Figure 3, purple shaded area) and disability weights. Internally consistent estimates of incidence and case fatality are derived from GBD data (by five-year age groups and sex) and Dismod II. Disability weights are derived from GBD disease specific YLDs and prevalence (by five-year age groups and sex) adjusted for all-cause YLDs (see section 2.1.1.4).
    Note that Dismod II is an external software. The code below is used to generate the inputs for Dismod II using GBD data (see Figure 3, purple shaded area). Dismod II requires data for the collection (population numbers and all cause mortality rates) and data set (disease specific parameters). We will generate a data set with all inputs for dismod by 5 year-age (represent at mid age) groups and sex for five diseases (ischemic heart disease, ischemic stroke, diabetes mellitus, colon and rectum cancer and breast cancer-women). WE WILL NEED ADJUSTMENT FOR TO MATCH RRS TO DISEASES (COLON AND RECTUM CANCER AND DIABETES MELLITUS). *Code here is work in progress**, for now the disease input data is organised externally with excel.

## Use externaly generated inputs for for disease life table. The previous code will serve to to some of the data preparation for DISMOD. Also need to devlop code to calculated DWs, for now, all generated in excel. 

idata <- read.csv("data/UK/idata.csv", stringsAsFactors = F)
## Use run_disease

i_disease <- c("ihd", "istroke", "diabetes", "colon_cancer", "breast_cancer")
i_sex <- c("male", "female")
i_age_cohort <- c(22, 27, 32, 37, 42, 47, 52, 57, 62, 67, 72, 77, 82, 87, 92, 97)

disease_life_table_list_bl <- list()
index <- 1

for (age in i_age_cohort){
  for (sex in i_sex){
    for (disease in i_disease) {
      # Exclude breast_cancer for Males
      if (sex == "male" && disease == "breast_cancer"){
       # cat("\n") #Uncomment to see list
      }
      else {
        # cat("age ", age, " sex ", sex, "and disease", disease, "\n") #Uncomment to see list
        disease_life_table_list_bl[[index]] <- run_disease(in_idata = idata, in_sex = sex, in_mid_age = age, in_disease = disease)
        index <- index + 1
      }
    }
  }
}
##### Uncommnet to check disease life table list
# View(disease_life_table_list_bl[[8]])

18. Generate mock change in incidence of disease to generate scenario life tables. In the final version, this will come from the calculated PIFs by disease. PIFs are applied here to incidence, however, may also be applied to case fatality, depending on the RRs. Also, we can include time delays from change in exposure to change in health outcomes via the PIF.

## Create value to use as factor changing incidence rates. 

incidence_change <- 0.95

#####Generate scenario incidence (for each disease)

incidence_sc <- list()
index <- 1

for (age in i_age_cohort){
  for (sex in i_sex){
    for (disease in i_disease) {
      
      # Exclude breast_cancer for Males
      if (sex == "male" && disease == "breast_cancer"){
        # cat("\n") # Uncomment to see list
      }
      else {
        
        incidence_sc[[index]] <- disease_life_table_list_bl[[index]]$incidence_disease * incidence_change
        index <- index + 1
        
      }
    }
  }
}

##### Uncommnet to check scenario incidence
# View(incidence_sc[[1]])

19. Use scenario incidence to calculate scenario disease life tables.

disease_life_table_list_sc <- list()
index <- 1
for (age in i_age_cohort){
  for (sex in i_sex){
    for (disease in i_disease) {
      # Exclude breast_cancer for Males
      if (sex == "male" && disease == "breast_cancer"){
        # cat("\n")
      }
      else {
        # cat("age ", age, " sex ", sex, "and disease", disease, "\n")
        # modify idata's incidence for the said scenario
        td1 <- idata
        td1[td1$age >= age & td1$sex == sex,][[paste("incidence", disease, sep = "_")]] <- incidence_sc[[index]]
        
        # Instead of idata, feed td to run scenarios
        disease_life_table_list_sc[[index]] <- run_disease(in_idata = td1, in_sex = sex, in_mid_age = age, in_disease = disease)
        disease_life_table_list_sc[[index]]$diff_inc_disease <- disease_life_table_list_sc[[index]]$incidence_disease - disease_life_table_list_bl[[index]]$incidence_disease
        disease_life_table_list_sc[[index]]$diff_prev_disease <- disease_life_table_list_sc[[index]]$px - disease_life_table_list_bl[[index]]$px
        disease_life_table_list_sc[[index]]$diff_mort_disease <- disease_life_table_list_sc[[index]]$mx - disease_life_table_list_bl[[index]]$mx
        disease_life_table_list_sc[[index]]$diff_pylds_disease <- (disease_life_table_list_sc[[index]]$px - disease_life_table_list_bl[[index]]$px) * disease_life_table_list_bl[[index]]$dw_disease
        index <- index + 1
      }
    }
  }
}
##### Uncommnet to check scenario life tables
# View(disease_life_table_list_sc[[1]])

20. Calculate life tables scenario.

#####Generate total change in mortality rate
##### Sum mortality rate scenarios (mx_sc_total)


mx_sc_total <- list()
l_index <- 1
index <- 1
for (age in i_age_cohort){
  for (sex in i_sex){
    mortality_sum <- NULL
    create_new <- T
    
    for (disease in i_disease) {
      if (sex == "male" && disease == "breast_cancer"){
        # cat("\n")
      }else{
        
        if (create_new){
          mortality_sum <- select(disease_life_table_list_sc[[index]], c('age', 'sex'))
          mortality_sum$total <- 0
          create_new <- F
          mortality_sum$total <- mortality_sum$total + (disease_life_table_list_sc[[index]]$diff_mort_disease)
        }else{
          mortality_sum$total <- mortality_sum$total + (disease_life_table_list_sc[[index]]$diff_mort_disease)
        }
        
        # cat(age, " - ", sex," - ",  disease," - ",  index, " - ", l_index,  "\n")
        index <- index + 1
      }
    }
    mx_sc_total[[l_index]] <- mortality_sum 
    l_index <- l_index + 1
  }
}  

##### Uncommnet to check sceanrio mortality and changes 
# View(mx_sc_total[[1]])


#####Generate total change in prevalent yld rates
#####total ylds rate= sum (change prevalence disease*dw)


pylds_sc_total <- list()
l_index <- 1
index <- 1
for (age in i_age_cohort){
  for (sex in i_sex){
    pylds_sum <- NULL
    create_new <- T
    
    for (disease in i_disease) {
      if (sex == "male" && disease == "breast_cancer"){
        # cat("\n")
      }else{
        
        if (create_new){
          pylds_sum <- select(disease_life_table_list_sc[[index]], c('age', 'sex'))
          pylds_sum$total <- 0
          create_new <- F
          pylds_sum$total <- pylds_sum$total +   (disease_life_table_list_sc[[index]]$diff_pylds_disease)
        }else{
          pylds_sum$total <- pylds_sum$total + (disease_life_table_list_sc[[index]]$diff_pylds_disease)
        }
        
        # cat(age, " - ", sex," - ",  disease," - ",  index, " - ", l_index,  "\n")
        index <- index + 1
      }
      
    }
    pylds_sc_total[[l_index]] <- pylds_sum
    l_index <- l_index + 1
  }
}  

##### Uncommnet to check scenario pyld change
# View(pylds_sc_total[[2]])

##### Calculate general life tables with modified mortality and pylds total#############################
###Original mortality rate is modified by the mx_sc_total (total change in mortality from diseases)
###Original pyld rate is modified by the change in each disease pylds


general_life_table_list_sc <- list()
index <- 1


for (age in i_age_cohort){
  for (sex in i_sex){
    # cat("age ", age, " and sex ", sex, "\n")
    # modify idata's mortality and pyld total for the said scenario
    td2 <- input_life_table
    td2[td2$age >= age & td2$sex == sex,][[paste("mx")]] <- general_life_table_list_bl[[index]]$mx + mx_sc_total[[index]]$total
    td2[td2$age >= age & td2$sex == sex,][[paste("pyld_rate")]] <- general_life_table_list_bl[[index]]$pyld_rate + pylds_sc_total[[index]]$total
    
    # Instead of idata, feed td to run scenarios
    general_life_table_list_sc[[index]] <- run_life_table(in_idata = td2, in_sex = sex, in_mid_age = age)
    
    index <- index + 1
  }
}
##### Uncommnet to check scenario life tables
# View(general_life_table_list_sc[[1]])
# View(general_life_table_list_bl[[2]])

21. Generate list of outputs by age and sex

#####In the following list "output_life_table", 32 data frames are nested per age and sex cohort

output_burden <- list()
l_index <- 1
index <- 1
for (age in i_age_cohort){
  for (sex in i_sex){
    
    #Males do not have breast cancer, that is why we need the if/else. 
    #We create a TRUE/FALSE variable for the loop to move into the next disease
    
    create_new <- T
    for (disease in i_disease) {
      if (sex == "male" && disease == "breast_cancer"){
        # cat("\n")
      }else{
        
        if (create_new){
          output_burden_sc <- select(disease_life_table_list_sc[[index]], c('age', 'sex', 'incidence_disease', 'mx', 'px'))
          names(output_burden_sc)[names(output_burden_sc) == 'incidence_disease'] <- paste('incidence_disease', disease, "sc", sep = "_")
          names(output_burden_sc)[names(output_burden_sc) == 'mx'] <- paste('mx', disease, "sc", sep = "_")
          names(output_burden_sc)[names(output_burden_sc) == 'px'] <- paste('px', disease, "sc", sep = "_")
          output_burden_bl <- select(disease_life_table_list_bl[[index]], c('incidence_disease', 'mx', 'px'))
          names(output_burden_bl)[names(output_burden_bl) == 'incidence_disease'] <- paste('incidence_disease', disease, "bl", sep = "_")
          names(output_burden_bl)[names(output_burden_bl) == 'mx'] <- paste('mx', disease, "bl", sep = "_")
          names(output_burden_bl)[names(output_burden_bl) == 'px'] <- paste('px', disease, "bl", sep = "_")
          
          ####New list to add calculations for changes in burden of disease (incidence and mortality numbers)
          
          output_burden_change <- list()
          
          output_burden_change$inc_num_bl <- disease_life_table_list_bl[[index]]$incidence_disease * (1 - disease_life_table_list_bl[[index]]$px) * general_life_table_list_bl[[l_index]]$Lx
          output_burden_change$inc_num_sc <- disease_life_table_list_sc[[index]]$incidence_disease * (1 - disease_life_table_list_sc[[index]]$px) * general_life_table_list_sc[[l_index]]$Lx
          output_burden_change$inc_num_diff <- (disease_life_table_list_sc[[index]]$incidence_disease * (1 - disease_life_table_list_sc[[index]]$px) * general_life_table_list_sc[[l_index]]$Lx) - (disease_life_table_list_bl[[index]]$incidence_disease * (1 - disease_life_table_list_bl[[index]]$px) * general_life_table_list_bl[[l_index]]$Lx)
          
          output_burden_change$mx_num_bl <- disease_life_table_list_bl[[index]]$mx * general_life_table_list_bl[[l_index]]$Lx
          output_burden_change$mx_num_sc <- disease_life_table_list_sc[[index]]$mx * general_life_table_list_sc[[l_index]]$Lx
          output_burden_change$mx_num_diff <- (disease_life_table_list_sc[[index]]$mx * general_life_table_list_sc[[l_index]]$Lx) - (disease_life_table_list_bl[[index]]$mx * general_life_table_list_bl[[l_index]]$Lx) 
          
          names(output_burden_change)[names(output_burden_change) == 'inc_num_bl'] <- paste('inc_num_bl', disease, sep = "_")
          names(output_burden_change)[names(output_burden_change) == 'inc_num_sc'] <- paste('inc_num_sc', disease, sep = "_")
          names(output_burden_change)[names(output_burden_change) == 'inc_num_diff'] <- paste('inc_num_diff', disease, sep = "_")
          names(output_burden_change)[names(output_burden_change) == 'mx_num_bl'] <- paste('mx_num_bl', disease, sep = "_")
          names(output_burden_change)[names(output_burden_change) == 'mx_num_sc'] <- paste('mx_num_sc', disease, sep = "_")
          names(output_burden_change)[names(output_burden_change) == 'mx_num_diff'] <- paste('mx_num_diff', disease, sep = "_")
          
          
          ###Bind all lists
          
          output_burden_sc <- cbind(output_burden_sc, output_burden_bl)
          output_burden_sc <- cbind(output_burden_sc, output_burden_change)
          
          create_new <- F
          
          #Here the calculations above are repeated, here is where the F is telling to move into the next disease
          
        }else{
          
          td3 <- select(disease_life_table_list_sc[[index]], c('incidence_disease', 'mx', 'px'))
          names(td3)[names(td3) == 'incidence_disease'] <- paste('incidence_disease', disease, "sc", sep = "_")
          names(td3)[names(td3) == 'mx'] <- paste('mx', disease, "sc", sep = "_")
          names(td3)[names(td3) == 'px'] <- paste('px', disease, "sc", sep = "_")
          
          td4 <- select(disease_life_table_list_bl[[index]], c('incidence_disease', 'mx', 'px'))
          names(td4)[names(td4) == 'incidence_disease'] <- paste('incidence_disease', disease, "bl", sep = "_")
          names(td4)[names(td4) == 'mx'] <- paste('mx', disease, "bl", sep = "_")
          names(td4)[names(td4) == 'px'] <- paste('px', disease, "bl", sep = "_")
          
          
          
          output_burden_change2 <- list()
          
          output_burden_change2$inc_num_bl <- disease_life_table_list_bl[[index]]$incidence_disease * (1 - disease_life_table_list_bl[[index]]$px) * general_life_table_list_bl[[l_index]]$Lx
          output_burden_change2$inc_num_sc <- disease_life_table_list_sc[[index]]$incidence_disease * (1 - disease_life_table_list_sc[[index]]$px) * general_life_table_list_sc[[l_index]]$Lx
          output_burden_change2$inc_num_diff <- (disease_life_table_list_sc[[index]]$incidence_disease * (1 - disease_life_table_list_sc[[index]]$px) * general_life_table_list_sc[[l_index]]$Lx) - (disease_life_table_list_bl[[index]]$incidence_disease * (1 - disease_life_table_list_bl[[index]]$px) * general_life_table_list_bl[[l_index]]$Lx)
          
          output_burden_change2$mx_num_bl <- disease_life_table_list_bl[[index]]$mx * general_life_table_list_bl[[l_index]]$Lx
          output_burden_change2$mx_num_sc <- disease_life_table_list_sc[[index]]$mx * general_life_table_list_sc[[l_index]]$Lx
          output_burden_change2$mx_num_diff <- (disease_life_table_list_sc[[index]]$mx * general_life_table_list_sc[[l_index]]$Lx) - (disease_life_table_list_bl[[index]]$mx * general_life_table_list_bl[[l_index]]$Lx) 
          
          names(output_burden_change2)[names(output_burden_change2) == 'inc_num_bl'] <- paste('inc_num_bl', disease, sep = "_")
          names(output_burden_change2)[names(output_burden_change2) == 'inc_num_sc'] <- paste('inc_num_sc', disease, sep = "_")
          names(output_burden_change2)[names(output_burden_change2) == 'inc_num_diff'] <- paste('inc_num_diff', disease, sep = "_")
          names(output_burden_change2)[names(output_burden_change2) == 'mx_num_bl'] <- paste('mx_num_bl', disease, sep = "_")
          names(output_burden_change2)[names(output_burden_change2) == 'mx_num_sc'] <- paste('mx_num_sc', disease, sep = "_")
          names(output_burden_change2)[names(output_burden_change2) == 'mx_num_diff'] <- paste('mx_num_diff', disease, sep = "_")
          
          
          ###Bind all lists
          
          output_burden_sc <- cbind(output_burden_sc, td3)
          output_burden_sc <- cbind(output_burden_sc, td4)
          output_burden_sc$age_cohort <- age 
          output_burden_sc <- cbind(output_burden_sc, output_burden_change2)
          
        }
        
        # cat(age, " - ", sex," - ",  disease," - ",  index, " - ", l_index,  "\n")
        index <- index + 1
      }
      
    }
    
    # general_life_table_list_sc and general_life_table_list_bl (Lx)
    output_burden_lf_sc <- select(general_life_table_list_sc[[l_index]], c('Lx', 'Lwx'))
    names(output_burden_lf_sc)[names(output_burden_lf_sc) == 'Lx'] <- paste('Lx', "sc", sep = "_")
    names(output_burden_lf_sc)[names(output_burden_lf_sc) == 'Lwx'] <- paste('Lwx', "sc", sep = "_")
    
    output_burden_lf_bl <- select(general_life_table_list_bl[[l_index]], c('Lx', 'Lwx'))
    names(output_burden_lf_bl)[names(output_burden_lf_bl) == 'Lx'] <- paste('Lx', "bl", sep = "_")
    names(output_burden_lf_bl)[names(output_burden_lf_bl) == 'Lwx'] <- paste('Lwx', "bl", sep = "_")
    
    
    output_burden_lf_sc$Lx_diff <- general_life_table_list_bl[[l_index]]$Lx - general_life_table_list_sc[[l_index]]$Lx
    output_burden_lf_sc$Lwx_diff <- general_life_table_list_bl[[l_index]]$Lwx - general_life_table_list_sc[[l_index]]$Lwx
    
    output_burden_sc <- cbind(output_burden_sc, output_burden_lf_sc)
    output_burden_sc <- cbind(output_burden_sc, output_burden_lf_bl)
    
    
    output_burden[[l_index]] <- output_burden_sc
    l_index <- l_index + 1
    
  }
}

#Uncomment to check
# View(output_burden[[32]])

22. Combine all lists of outputs in a data frame to facilitate extraction of outcomes of interest and plotting.

#####Generate a data frame for all results and create function to get outcomes. 

output_df <- plyr::ldply(output_burden, rbind)

#Remove variables that are not used in the generation of outputs. CHANGE THIS NAMES, TOO LONG

output_df <- subset(output_df, select = -c(incidence_disease_ihd_bl, incidence_disease_ihd_sc, incidence_disease_istroke_bl, incidence_disease_istroke_sc, incidence_disease_diabetes_bl, incidence_disease_diabetes_sc, incidence_disease_breast_cancer_bl, incidence_disease_breast_cancer_sc, 
incidence_disease_colon_cancer_bl, incidence_disease_colon_cancer_sc,  mx_ihd_bl, mx_ihd_sc, mx_istroke_bl, mx_istroke_sc, mx_diabetes_bl, mx_diabetes_sc, mx_breast_cancer_bl, mx_breast_cancer_sc, mx_colon_cancer_bl, mx_colon_cancer_sc, px_ihd_bl, px_ihd_sc, px_istroke_bl, px_istroke_sc, px_diabetes_bl, px_diabetes_sc, px_breast_cancer_bl, px_breast_cancer_sc, px_colon_cancer_bl, px_colon_cancer_sc))

23. Plot outcomes for each age and sex cohort using created function. ADD A LOOP TO GENERATE ALL OUTCOMES FOR IN_AGE, IN_POPULATION, IN_OUTCOMES AND IN_LEGEND. FIX THE COLORS IN THE FUNCTION.

## Do a loop, the first in_outcomes is age.

i_outcome <- c("mx", "inc")

for (age in i_age_cohort){
  for (sex in i_sex) {
     for (disease in i_disease){
       for (outcome in i_outcome) {
         {
       if (sex == "male" && disease == "breast_cancer"){
         # cat("\n")
       }else{

 p <- plot_output(in_data = output_df, in_age = age, in_population = sex, in_outcomes = c("age", paste(outcome, "num", "bl", disease, sep = "_"), paste(outcome, "num", "sc", disease, sep = "_"), paste(outcome, "num", "diff", disease, sep = "_")), in_legend = paste(outcome, disease, sep = " "))
           }
         }
       }
     }
   }
 }

paste(age, sex, disease, outcome, sep = “_“)

24. Add up all outcomes per year of simulation with generated function.

####Generate data frame with all outputs to graph total change in burden by simulation year. 
###first, need to run function for males and females separetly, in_cohorts indicates the number of age cohorts to 
####include. To show all select 16. what is specified in in_outcomes will be graphed or can also be presented as a total. 
aggregate_frame_males <- gen_aggregate(in_data = output_df, in_cohorts = 16, in_population = "male", in_outcomes = c('inc_num_bl_ihd','inc_num_sc_ihd', "inc_num_diff_ihd"))

aggregate_frame_females <- gen_aggregate(in_data = output_df, in_cohorts = 16, in_population = "female", in_outcomes = c('inc_num_bl_ihd','inc_num_sc_ihd', "inc_num_diff_ihd"))

#####The following adds up both males and females
# Remove non-numeric columns starting with age and sex
aggregate_frame_males <- aggregate_frame_males %>% select(-starts_with("age"), -starts_with("sex"))

aggregate_frame_females <- aggregate_frame_females %>% select(-starts_with("age"), -starts_with("sex"))

# Create a copy of aggregate_frame_females
total_aggr <- aggregate_frame_females
# Add aggregate_frame_males values to it
for (i in 1:ncol(aggregate_frame_females)){
  total_aggr[i] <- total_aggr[i] + aggregate_frame_males[i]
}
total_aggr$sim_year <- seq.int(nrow(total_aggr))

25. Plots using total aggregated change in burden.

####This plot has to be customised to in_outcomes, here, only totals shown, but specifications are up to the user

####[] is used here to indicate the number of simulation years into the future. 
####Disease outcomes has to be changed to the outcome of interest

disease_outcome <- "Deaths ischemic heart disease"

total_plot <- ggplot(total_aggr[1:20,], aes(x = sim_year)) +
  
  geom_line(mapping = aes(y = total_inc_num_bl_ihd, colour = "total_mx_num_bl_ihd")) +
  theme_classic() +
  geom_hline(yintercept=0, linetype="dashed", color = "black") +
  geom_line(mapping = aes(y = total_inc_num_sc_ihd, colour = "total_mx_num_sc_ihd")) +
  geom_line(mapping = aes(y = total_inc_num_diff_ihd, colour = "total_mx_num_diff_ihd")) +
  xlab ("Simulation years") + ylab ("Cases") + labs (title = paste(disease_outcome)) +
  theme(plot.title = element_text(hjust = 0.5, size = 12)) +
  scale_color_discrete(name = paste(""), labels = c("Baseline", "Difference", "Scenario")) +
  theme(plot.title = element_text(hjust = 0.5))
  

#Print to view
print(total_plot)